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Abstract 

The universal scaling function of the square lattice Ising model in a magnetic field is 
obtained numerically via Baxter's variational corner transfer matrix approach. The high 
precision numerical data is in perfect agreement with the remarkable field theory results 
obtained by Fonseca and Zamolodchikov, as well as with many previously known exact and 
numerical results for the 2D Ising model. This includes excellent agreement with analytic 
results for the magnetic susceptibility obtained by Orrick, Nickel, Guttmann and Perk. In 
general the high precision of the numerical results underlines the potential and full power of 
the variational corner transfer matrix approach. 
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The Ising model has played a prominent role in the development of the theory of phase 
transition and critical phenomena [1-8]. The partition function of the nearest-neighbour Ising 
model on the square lattice reads 

Z = ^exp^P^aiaj + H^ai^ , (jj = ±1, (1) 

where the first sum in the exponent is taken over all edges, the second over all sites and the outer 
sum over all spin configurations {u} of the lattice. The constants H and (3 denote the (suitably 
normalized) magnetic field and inverse temperature. The specific free energy, magnetization and 
magnetic susceptibility are defined as 

1 dF d'^F 

where is the number of lattice sites. The model exhibits a second order phase transition at 
P = H = 0, where [1] 

= i log(l + V2) = 0.44068679 .... (3) 

In what follows we will exclude the temperature variable /? in favour of a new variable 

2t = cosech(2/?) - sinh(2/3) , = 0, (4) 

which is vanishing for (3 = Pc and positive for P < Pc (above the critical temperature). Note also 
that this variable changes sign under the Kramers- Wannier duality transformation for H = 0. 
Another useful related variable is 



k = k{T) = (Vl + r2 + r)2. (5) 

According to the scaling theory [6,9,10], the leading singular part, Fsing{T, H), of the free 
energy ([2]) in the vicinity of the critical point can be expressed through a universal function 
J^{m, h), 

Fsing{T,H)=r{m{T,H),h{T,H)), T ^ 0, H ^ 0, (6) 

where r and H enter the RHS only through non-linear scaling variables [11], 

m = m{T, H) = -V2t + ©(r^) + 0{H^) + . . . , h = h{T, H) = ChH + HO{T) + 0{H^) . . . , 

(7) 

which are analytic functions of r and H. The coefficients in these expansions (for instance, the 
leading coefficients — \/2 and Ch) are specific to the square lattice Ising model, however, the 
function T{m, h) is the same for all models in the 2D Ising model universality class. It can be 
written as 

.F(m,/i) = ^logm2 + /ii6/i5$(^)^ r^=^ (8) 

where ^{rf} is a universal scaling function of a single variable r/ (the scaling parameter). 

The function T{m, h) has a concise interpretation in terms of 2D Euclidean quantum field 
theory. Namely, it coincides with the vacuum energy density of the "Ising Field Theory" (IFT) 
[12]. The latter is defined as a model of perturbed conformal field theory with the action 

^IFT = -4(c=i/2) + ^ y ^(^) d'^x + h j (^i^) d'^x , (9) 
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where .4.(^=1/2) stands for the action of the c = 1/2 CFT of free massless Majorana fermions, 
(t{x) and e{x) are primary fields of conformal dimensions 1/16 and 1/2. Their normahzation is 
fixed by the usual CFT convention, 



x|2(e(x)e(0)) ^1; 



X 



{a{x)a{0)) ^ 1 as |a;|^0. 



(10) 



With this normalization, the parameters m and h have the mass dimensions 1 and 15/8, respec- 
tively, and the scaling parameter 77 in (l8|) is dimensionless. 

The scaling function ([8]) is of much interest as it controls all thermodynamic properties of 
the Ising model in the critical domain. Although there are many exact results (obtained through 
exact solutions of ([9]) at /i = and all r [1,3,13-16], and at r = and all h [8,17-22]; these data 
are collected in [23]) as well as much numerical data [24-28] about this function, its complete 
analytic characterization is still lacking. 

Recently [12] the function particularly its analytic properties, have been thoroughly 
studied in the framework of the Ising Field Theory ([9]). The authors of [12] made extensive 
numerical calculations of the scaling function $(r/) using the "Truncated Free-Fermion Space 
Approach" (TFFSA), which is a modification of the well-known "Truncated Conformal Space 
Approach" (TCSA) [29,30]. 

The primary motivation of our work was to confirm and extend the field theory results 
of ref. [12] through ab initio calculations, directly from the original lattice formulation ([T|) of 
the Ising model. We used Baxter's variational approach based on the corner transfer matrix 
method [31-33]. The main advantage of this approach over other numerical schemes (e.g., 
the row-to-row transfer matrix method) is that it is formulated directly in the limit of an 
infinite lattice. Its accuracy depends on the magnitude of truncated eigenvalues of the corner 
transfer matrix (which is at our control), rather than the size of the lattice. The details of our 
calculations along with numerical data for the free energy, magnetization and internal energy 
of the Ising model will be presented elsewhere [34]. We used several important enhancements 
of the original Baxter approach [32], in particular an improved iteration scheme [35], known as 
the corner transfer matrix renormalization group (CTMRG). The calculations were performed 
for a grid of values of the magnetic field and temperature in the range 10~^ < H < 10~^ and 
0.9/3c < /? < ^-^fic, containing about 10,000 distinct point (excluding a small region around the 
critical point). The results for the scaling function ^{rj) are shown in Fig. [TJ All calculated 
data points collapse on a smooth curve, shown by the solid line (additional details presented on 
the picture are explained below). Our numerical results for $(r/) remarkably confirm the field 
theory calculations of [12], to within all six significant digits presented thereirH. 

The precision of our numerical calculations was tested against all available exact results for 
the Ising model. In particular, the agreement between calculated and exact values for the zero- 
field free energy [1], magnetization [2] and magnetic susceptibility [36] in our working range of 
temperatures varied between 14 and 28 decimal places (depending on the distance to the critical 
point). In addition to these checks we also confirmed and extended many previously existing 
numerical results for the Ising model. Some details of our results are described below. 

For the following discussion it is convenient to rewrite ([8]) in an alternative form, introducing 
two more scaling functions. 



^We thank Alexander Zamolodchikov for providing us with additional unpublished numerical data for ^{rj), 
which are again in perfect agreement with our results. 




77i < 
m > 




(11) 
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where Ghigh{Q) = Giow{0) = 0, corresponding to 

2 

jr(m,0) = — logm^ . (12) 

These scaling functions are thoroughly discussed in [12]. The function GhighiO can be expanded 
in a series in even powers of ^ 

GMghiO = G2f + G4^^ + Gef + ... (13) 

convergent in some domain around the origin of the ^-plane. The function (jr;ou)(0 admits an 
asymptotic expansion 

Giov^iO = GiC + + Gse + ■■■ (14) 

for small positive ^. These new functions are simply related to $(?7). Note, in particular, that 
the coefficients G„ and G„ control the behavior of the function $(77) for large values of r) on the 
real line, 

1 — 7 — 29 

^lowii]) = Giijs + G2r]~i + G3r]~~a + ... for real rj +00 , (15) 

'^highiv) = G2{-vy^ +G4i-vy^ +G6{-r]y^ + ... for real r? ^ -00 . (16) 
Finally, for small values of rj, 

2 °° 

Hv) = -^^og7j' + Y,'^kv' (17) 

fe=0 

where the series converges in a finite domain around the origin of the complex 77-plane. 

Some of the above expansion coefficients are known exactly. The coefficient Gi is known 
explicitly [4] 

Gi = -2^/12 j3/2 ^ _i.357838341706595 . . . , (18) 

where A = 1.282427... is the Glaisher constant. The coefficients G2 and G2 have integral 
expressions [14, 15] involving solutions of the Painleve III equation. They were numerically 
evaluated to very high precision (50 digits) in [36], 

G2 = -1.845228078232838 . . . , 62 = -0.0489532897203 .... (19) 

The coefficient $0 was calculated in [17], 

^ _ (27r)T^7(5)7(|)7(n) ^ _i.i9773338379799 . . . , ^(x) = (20) 

where T{x) is the standard F-function. The coefficient $1 has an explicit integral representation, 
obtained in [37]. We have evaluated the required integral explicitly, 

32 2! ^(^h(|) n 7(^) 

$1 = 7 ^"^ 19 = -0.3188101248906 .... (21) 

225 (27r)T5 [^(i)^2(^)]T5 

Contrary to the field theory case, the lattice free energy, 

F{t,H) H) + Freg{T, H) + Fsub{T, H), T,H^O, (22) 
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never coincides with its leading universal part It also contains regular terms Freg{T, H), 
analytic in r and H, as well as subleading singular terms Fsui,{t, H), which are non-analytic, but 
less singular than the first term in (122^ . Therefore, to extract the universal scaling function from 
the lattice calculations one should be able to isolate and subtract these extra terms. Moreover, 
one needs to know the explicit form of the non- linear scaling variables ([7]). In principle, all 
this information can be determined entirely from numerical calculations (provided one assumes 
the values of exponents of the subleading terms, predicted by the analysis [36, 38] of the CFT 
irrelevant operators, contributing to the free energy (I22p ). The accuracy of this "fully numerical" 
approach, however, deteriorates rapidly for the higher order terms. Much more accurate results 
can be obtained if the numerical work is combined with known exact results. 
Write the non-linear variables ([7]) in the form. 



c(t) + d{T) + 0{H^ 



H) = -V2t a{T) + H^b{T) + 0{H^ 
h{T,H) = ChH 

where a(0) = c(0) = 1, /i(t, H) = —h{T, —H). Similarly, write the regular part as, 

Freg{T, H) = A{t) + B{t) + 0{H^) . 



(23) 



(24) 



As shown in [36], the most singular subleading term, contributing to (I22p is of the order of 
0(t^/^ H^) (see Eq. (|36p below). The fact that subleading terms arise only in such a high order 
(~ m^) is a remarkable property of the square lattice nearest-neighbour Ising model, which 
greatly simplifies the calculation of the universal scaling function from the numerical data. 
For H = the expression ([22]) should reduce to Onsager's exact result [1] 



F{t, 0) = log ^2 cosh(2/3) + / — log 



1+1 



cos 



1 + 



1/2- 



(25) 



Therefore one should be able to rewrite the last formula in the form (I12p plus regular terms. 
This is achieved by choosing [39] 



air) 



/'dxF(i,i,l;-xr2)/(l+xr2)V2 
Jo 



1/2 



3r2 137r^ 
le" 1536 



+ 0(t^ 



(26) 



where F{a, b, c, z) denotes the Gauss hypergeometric function. The corresponding contribution 
to (EH) then reads 



Air) = 

vr 



log 2 T r2(l-K51og2) 5r^(l 6 log 2) 



+ 2 



12 + 



647r 



+ 0(tS 



(27) 



where Q = 0.915965594 ... is the Catalan constant. 

Next, using the exact expression for the zero field magnetization [2,3] 

M(t,0) = (1 - A:(r)2)i/8^ r < 
with /c(r) defined in ([5]), one finds from (P, ([U]), ^ and ([22]) 

, . M(r) _ r 15r2 9r^ 4333r'^ 



c r 



(-4ra(r))V8 



, r 15r^ 
^+4 + 12^ 



512 98304 



+ Oir^ 



(28) 



(29) 



5 



and also 

Ch = -23/^V<5i = 0.838677624411 .... (30) 

Finally, consider the zero-field susceptibility. No simple closed form expression for the zero- 
field susceptibility x{t) is known. However, the authors of [36] obtained remarkable asymptotic 
expansions of x(t) for small r to within 0(r^^) terms with high-precision numerical coefficients. 
We have compared their series with our numerical results for the susceptibility in the range of 
temperatures |r| = 0.05-0.14 and found that they agree to each other in 14 to 18 significant 
digits (depending on the value of r). Their result can be written as (retaining the terms up to 
0(|r|^/^), inclusive) 

x(r)oNGP = -2-^^^ G"(0) 1-1-^1 + i + ^ + + 

+e(T) + fir) log |r| + 0{t^ log |r|) (31) 

where e(r) and /(r) are explicitly known second order polynomials in r. The coefficient in front 
of the ffist term is written in our notation^. The symbol G"(0) there stands for the second 
derivative of the scaling function G(^) with respect to its argument. Namely G"{0) = G'^igfiiO) = 
2G2 for r > and G"(0) = ^^'^^^(O) = 2G2 for r < 0. The corresponding values are given in 
(fT9]) above. Next, calculating the second field derivative of (f22]l at -fT = 0, one obtains 



Xir) = -2-iClG"{0)\T\~ia{T)-UiT 



H=0 



-25(r) + :^^i^(l + log(2rMr))) • (32) 



vr 



Equating this expression to (j3ip and using (j26p . (j29p and explicit forms of the polynomials e(r) 
and /(t) from [36], one obtains 

B{t) = 0.0520666225469 + 0.0769120341893 r + 0.0360200462309 + 0{t^) , (33) 

and 



Noting that 



6(r) = i_ih(l + ^ + 0(r2)) , fih = 0.071868670814 . (34) 

air) 4 c(tY = H \ \ h O(t^) (35) 

\ J \ J 2 8 16 192 ^ ^ ^ ^ 

one immediately obtains the main contribution to the subleading term 

2-ic^G"(0))-'5!%5ijiQ „ ^ = -4Mf+... . (36) 



H=o 384 



The results of [36] provided the first convincing demonstration of the violation to simple one- 
parametric scaling in the square lattice Ising model. Note that despite the 0(|t|^/^) term gives 



^The correspondence with the notations of ref. [36] is as follows. Their Pxij) is denoted here as x(''')ongp- 
Their coefficients C± are connected to our constants by 

C+ = -2^ (2/3eV2)"5 CIG2 , C- = -2^ (2/3,72)^5 , 

where l3c is given in ([3|. 
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a very small contribution to the susceptibility we were able to confidently quantify it from our 
numerical results. Namely, we estimated the coefficient of the term in the series shown in 
parenthesis in the first line of (j3ip . Our estimate is 23.004/384 which is within 0.02% of its 
exact value 23/384. 

The coefficient (i(r) in (123p was estimated from our numerical data for the internal energy, 

d{T)=eh + 0{T), eft = -0.007(1), (37) 

which is in agreement with the result Ch = —0.00727(15) from [26]. 

The above expressions were used to analyze our extensive numerical data and extract the 
necessary information to obtain the universal scaling function. The results are summarized 
in four tables. For convenience of comparison we quoted the corresponding results from [12], 
including those obtained through the TFFSA, the high/low-temperature and extended disper- 
sion relations (DR). Earlier exact and numerical results for the same quantities are also quoted 
(whenever available). Figure 1 shows about 10,000 data points for the scaling function $(r/). As 
expected, the points collapse on a smooth curve (their spread is much smaller than the resolution 
of the picture). Figure 1 also shows plots of the asymptotic expansions p3]) . and (fT7|) with 
the coefficients given in Table 1, Table 2 and Table 3. These expansions are seen to "stitch" 
together very well and give a reasonably good analytic approximation to ^{rj) in the whole real 
line of r/. Our Fig. [T] essentially coincides with Fig. 10 of [12]. The numerical values of ^{rj) at 
small integer values of r] are given in Table 4. 

The calculations were performed on the 24-processor Linux Cluster System at the ANU 
Research School of Physics and Engineering and on the 1928-processor SGI Altix 3700 Bx2 
Cluster at the ANU Supercomputer Facility. The level of parallelization varied between 15 and 
50 processors. The total amount of CPU time spent for this work was about 9,000 hours (single 
processor equivalent). 

To conclude, we have implemented Baxter's variational corner transfer matrix approach to 
obtain the universal scaling function of the square lattice Ising model in a magnetic field, as 
shown in Figure 1 and Table 4. The numerical data is seen to be in remarkable agreement with 
the field theory results obtained by Fonseca and Zamolodchikov [12]. We also report a remarkable 
agreement (11 to 14 digits) between our numerical values for Gi, G2 and G2 and the classic 
exact results of Barouch, McCoy, Tracy and Wu [4,14,15] (see Tables 1 and 2), and a similar 
agreement between the values $0 and $1 and the exact predictions [17,37] of Zamolodchikov's 
integrable field theory [8] (see Table 3) . The high precision of the numerical results underline 
the full power and further potential of the variational corner transfer matrix approach. In this 
case, the results show beyond any doubt the validity of the connection between the scaling limit 
of the Ising model in a magnetic field and the Ising Field Theory ([9]). 

The authors thank H. Au-Yang, R.J. Baxter, G. Delfino, S.L. Lukyanov, J.H.H. Perk, 
C.A. Tracy and A.B. Zamolodchikov for useful discussions and remarks. This work has been 
partially supported by the Australian Research Council. 
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Figure 1: Scaling function of the two-dimensional Ising model in a magnetic field. The figure shows 
about 10,000 data points for the scaling function $(r/). As expected, the points collapse 
on a smooth curve (their spread is much smaller than the resolution of the picture). 
The figure also shows plots of the asymptotic expansions ([15]), ((T6)) and (fTT)) with the 
coefficients given in Table 1, Table 2 and Table 3. 
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CTM (This work) 


High-T DR [12] 


From References 


G2 


-1.8452280782328(2) 


-1.8452283 


-1.845228078232838... 


[14,15,36] 


G4 


8.333711750(5) 


8.33410 


8.33370(1) 


[26] 


Ge 


-95.16896(1) 


-95.1884 


-95.1689(4) 


[26] 


Gs 


1457.62(3) 


1458.21 


1457.55(11) 


[26] 




-25891(2) 


-25889 


-25884(13) 


[26] 



Table 1: Numerical values of the coefficients G2n in (I13p . The second column contains the high- 
temperature dispersion relation (DR) results from [12]. The coefficient G2 is known exactly 
[14,15,36]. 





CTM (This work) 


Low-T DR [12] 


From References 




-1.3578383417066(1) 


-1.35783835 


-1.357838341706595... [4] 


G2 


-0.048953289720(1) 


-0.0489589 


-0.0489532897203... [14,15,36] 


G3 


0.038863932(3) 


0.0388954 


0.0387529 [40]; 0.039(1) [25] 


G4 


-0.068362119(2) 


-0.0685060 


-0.0685535 [40]; -0.0685(2) [25] 


G5 


0.18388370(1) 


0.18453 




Ge 


-0.6591714(1) 


-0.66215 




G7 


2.937665(3) 


2.952 




Gs 


-15.61(1) 


-15.69 





Table 2: Numerical values of the coefficients in (fT4| . The second column contains the low- 
temperature dispersion relation results from [12]. The right column refers to exact values 
of Gi [4] and G2 [14,15,36] and other numerical results. 
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CTM (This work) 


TFFSA [12] 


Ext. DR [12] 


From References 




-1.197733383797993(1) 


-1.1977331 


-1.1977320 


-1.19773338379799339... [17] 


$1 


-0.318810124891(1) 


-0.3188103 


-0.3188192 


-0.3188101248906... [37] 


$2 


0.110886196683(2) 


0.1108867 


0.1108915 




$3 


0.01642689465(2) 


0.0164266 


0.0164252 




$4 


-2.639978(1) x 10"^ 


-2.64 X 10-^ 


-2.64 X 10-"^ 




$5 


-5.140526(1) X 10"^ 


-5.14 X 10-^ 


-5.14 X 10-"^ 




$6 


2.08865(1) X 10"^ 


2.07 X 10""^ 


2.09 X 10"^ 




<I>7 


-4.4819(1) X 10-5 


-4.52 X 10-5 


-4.48 X 10-5 




^>8 






3.16 X lO-'^ 




$9 






4.31 X 10-^ 




^10 






-1.99 X 10"^ 





Table 3: Numerical values of the coefficients $„ in (fT7|) . The second and third columns contain 
results from [12], obtained through the TFFSA and the extended dispersion relations 
(DR), respectively. The forth column refers to exact results; the numerical value of <I>i 
therein is taken from Eq. (PT|) . 



$(r/) 


CTM (This work) 


TFFSA [12] 


High/Low-T DR [12] 


Ext 


DR [12] 


$(-5) 


-0.10920919 


-0.1092101 


-0.1092092 


-0 


1088626 


$(-4) 


-0.15926438 


-0.1592682 


-0.1592643 


-0 


1589421 


$(-3) 


-0.25298908 


-0.2529928 


-0.2529887 


-0 


2527417 


$(-2) 


-0.44132564 


-0.4413450 


-0.4413249 


-0 


4412136 


$(-1) 


-0.78396650 


-0.7839665 


-0.7839668 


-0 


7839576 




-1.19773338 


-1.1977330 




-1 


1977320 




-1.38984135 


-1.3898410 


-1.3898417 


-1 


3898063 


m 


-1.49305602 


-1.4930558 


-1.4930566 


-1 


4929849 


m 


-1.56427320 


-1.5642732 


-1.5642736 


-1 


5641727 


m 


-1.61885066 


-1.6188506 


-1.6188510 


-1 


6187275 


$(5) 


-1.66324828 


-1.6632483 


-1.6632485 


-1 


6631076 



Table 4: Numerical values of <&(?7) at small integer values of 77. The second, third and fourth 
columns contain results [12] from the TFFSA, the high/low-temperature dispersion and 
the extended dispersion relations (DR) . 
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